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THE TURBULENT HEAT FLUX IN LOW MACH NUMBER FLOWS 


WITH LARGE DENSITY VARIATIONS 


Peter J. O’Rourke and Lance R. Collins 
Theoretical Division, Group T-3 
Los Alamos National Laboratory 
Los Alamos, New Mexico 87545 


I. INTRODUCTION: THE DIRECTED ENERGY FLUX 

This paper is concerned with a physical effect of fundamental importance 
for the modeling of turbulence transport in flows with large density variations. 
The effect occurs because the interaction of pressure and density gradients gives 
rise to a turbulent heat flux, which we call a directed flux, that is not accounted 
for in turbulence models for constant density flows. To see how this flux arises, it 
is perhaps best to consider an example of Rayleigh-Taylor instability, as depicted 
in Fig. 1. A heavy, cold gas overlays a light, hot gas in a box with a gravitational 
acceleration in the negative z-direction. The induced hydrostatic pressure gradi- 
ent accelerates the light gas into the heavy gas and causes the instability and 
mixing. The velocities averaged across a horizontal plane are in the z-direction, 
and because the heavy gases are falling, the mass-averaged velocities in the 
mixing region will be negative. Relative to a surface moving downward with the 
mass-averaged velocity, there will be a net upward flux of energy. This is 
because although the mass flux of light gas crossing the surface upward equals 
the mass flux of heavy gas crossing downward, the light gas, being hotter, carries 
with it more energy per unit mass. 

This upward energy flux is the directed flux. In the example of Fig. 1, this 
energy flux is in the direction of the negative of the temperature gradient, just as 
given by the laminar Fourier heat conduction law. If after some time we were to 
turn the box over so that the light gas overlay the heavy gas, to the extent that 
the two gases had not already mixed on the molecular level, there would be an 
unmixing in which the light gas would separate from the heavy gas. In this 
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Fig. 1. Schematic depiction of 
Rayleigh-Taylor instability. 
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unmixing process, the directed heat flux would be upward, in the direction of the 
mean temperature gradient and opposite the direction given by the Fourier heat 
conduction law. This phenomenon has thus been called countergradient diffu- 
sion. l It cannot be predicted by turbulence models that use gradient transport, or 
a Fourier-like law, to describe turbulent heat transport. In our example, the tur- 
bulent heat flux <t>^ is in the direction opposed to the pressure gradient, rather 
than the temperature gradient. We shall see that taking <j>/i ~ — vp is often more 
realistic for turbulent flows. 

Between single-phase, two-density turbulent flows and two-phase flows 
there is an analogy that we will exploit in our turbulence modeling. This analogy 
will be used in helping to formulate the equations and in developing closure 
approximations for some of the terms. In two-phase flow modeling, which has 
received much attention within the last ten years 2 . 3 separate mass, momentum, 
and energy equations are kept for each phase, and these equations are coupled 
through functions that give the exchange rates of mass, momentum, and energy 
between the phases. Following Besnard, Harlow, and Rauenzahn, 4 we will use 
an alternative, analogous formulation. In place of two mass equations, we will 
keep an equation for the mean density and one for density fluctuations. In place 
of two momentum equations, we will keep a mean momentum equation and an 
equation for mean velocity differences associated with fluid elements of differing 
density. We use this second formulation because it allows for the possibility of 
modeling not just two-density flows, but the flows with a spectrum of densities 
that often occur in practical applications. 

Two physical examples, one of two-phase flow and one of single-phase flow 
with density variations, serve to illustrate the analogy and another situation in 
which countergradient transport can arise. In both examples, pressure gradients 
are responsible for centrifuging lighter material inward toward the centers of 
rotating flows. In a two-phase bubbly flow, this effect has been observed in the 
vortices in the wake of an obstacle. 3 In single-phase flow, it is probably responsi- 
ble for the experimental results of Wahiduzzaman and Ferguson.6 The experi- 
menters measured the radial temperature profiles in an axisymmetric swirling 
flow in a constant volume cylinder. The experimentally measured temperatures 
are plotted at four different times as the circles in Fig. 2. The lines are computed 
temperature profiles using the KIVA code? with a k-z turbulence model® and 
gradient heat transport with a turbulent Prandtl number of 0.9. It can be seen 
that a hot region in the center of the cylinder persists much longer in the experi- 
ment than in the calculation, showing the large errors that can arise when a 
gradient heat transport approximation is used. 

The phenomenon of countergradient transport in single-phase flows was 
recognized seven years ago in research on the structure of turbulent premixed 
flames. 1 In retrospect, it is easy to see how this phenomenon arises. Figure 3 
depicts schematically a planar turbulent premixed flame with velocities shown in 
the frame of reference of the flame. Mass conservation and the fact that the com- 
bustion is nearly isobaric together imply that the hot product gas velocities will 
be larger than those in the reactants. Momentum conservation then implies that 
the pressure in the products will be lower than in the reactants. Since the dir- 
ected heat flux is in the direction opposed to the pressure gradient, this heat flux 
will be from the colder to the hotter gases; that is, it will be countergradient 
transport. 

Two approaches have been used for modeling turbulent premixed flames — a 
single-phase formulation and a two-phase formulation. In the single-phase for- 
mulation of Bray, Moss, and Libby 1 . 311 (BML formulation), equations are kept 
for the mean product gas concentration, the mean momentum, the turbulent 
fluxes of these quantities, and for the dissipation rate of turbulent kinetic energy. 
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Fig. 3. Schematic depiction of a turbulent premixed flame. 

Spaldingl2 utilizes a two-phase formulation, retaining mass, momentum, and 
energy equations for each phase. Spaldingl2 assumes that the major source of 
mixing is due to the difference between the mean velocity of the phases and thus 
ignores the turbulent kinetic energy within each phase. The BML formulation 
accounts for both sources to the turbulent kinetic energy. In principle, the equa- 
tions of one formulation should be derivable in terms of those of the other, al- 
though to our knowledge such a derivation and comparison have not been made. 
In this paper some closure approximations are proposed, based on a derivation of 
the single-phase equations from two-phase equations. Only the BML formulation 
has been compared with experimental measurements of turbulent flames, l 1 and 
satisfactory agreement was obtained. 

In practical applications of turbulent combustion, other physical effects that 
cause mixing and unmixing are superimposed on the pressure drop across the 
flame. At Los Alamos, we have been involved for the past twelve years in the 
numerical modeling of combustion in internal combustion engines. 7,13-16 
Figure 4 illustrates some of the complexities of the turbulence/chemistry inter- 
action in an engine burning premixed charges. A turbulent premixed flame is 
propagating away from an ignitor located near the center of the cylinder head 
wall. Mach numbers are small, and thus the mean pressure is nearly uniform in 
spaced and changing with time due to piston motion, combustion, and wall heat 
loss. Near the top of its motion, the piston, and the axial flow velocities in the 
combustion chamber, decelerate. This causes a small positive axial pressure 
gradient and induces Rayleigh-Taylor instability and mixing where the flame is 
propagating downward in the axial direction. This same pressure gradient will 
cause a differential axial acceleration of the hot products and cool reactants and 
promote Kelvin-Helmholtz instability where the flame is propagating radially. 
Swirl, a nearly symmetric rotational motion of the burning gases, is introduced 
by engine designers to promote mixing but will have two competing effects in the 
engine of Fig. 4. Swirl induced shears will enhance turbulence and mixing, but 
the radial pressure gradient caused by the swirl will, as in the experiments of 
Wahiduzzaman and Ferguson, 6 cause countergradient transport and suppress 
mixing. It is important to point out that among these various turbulence effects, 
only those associated with shear instability and mixing are accounted for in 
current engine models. 

In predicting turbulence in internal combustion engines and other practical 
combustors, one cannot use two-phase models or single-phase equations for two 
density flows. This is because within the reactants and products there will be 
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Fig. 4. The turbulence/chemistry interaction in an internal combustion engine. 

distributions of density. In the engine of Fig. 4 density differences will arise due 
to wall heat transfer and due to entropy differences in the product gases of this 
confined burn. 18 In other combustors density distributions are caused by charge 
non-uniformity and spray vaporization. 

In the next section we derive preliminary equations for a turbulent fluid 
with large density variations. Our aim is to develop a model that has three 
attributes: 

(1) the model can predict mixing and unmixing due to shear instabilities and 
pressure-density gradient interactions; 

(2) the model can account for a distribution of densities; and 

(3) the model equations can be efficiently integrated in two and three 
dimensions. 

The second attribute precludes use of two-phase flow equations, although investi- 
gation of the two-density limit will yield valuable information. The third attri- 
bute precludes use of the Reynolds stress equations, especially in three dimen- 
sions. It seems appropriate to seek a simple one- or two-equation extension of 
popular two-equation models for turbulent shear flows. 

II. THE EQUATIONS 

A. Overview 

We first derive equations for the average density p and the Favre-averaged 
velocity u and enthalpy h. Using the low Mach number assumption, we relate 
the turbulent heat flux to the difference between the average velocity u and 
the Favre-averaged velocity u. We denote this difference u — u by a, and the 
transport equation for a is derived and discussed. Closure approximations for 
terms in the a-equation are postulated based on the analogy between two-phase 
flows and single-phase, two-density flows. A comparison between the single- 
phase and two-phase equations suggests that the fluctuating stress terms in the 
a-equation are primarily associated with the decay of a. We present an algebraic 
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closure approximation for a that results in a heat flux that is the sum of contribu- 
tions due to gradient and directed transport. Our a-equation is compared with 
similar equations in the literature. 

B. The Equations of a Low Mach Number Flow with Large Density Variations 

For the low Mach number flow of a single component ideal gas with large 
density variations, the equations are the following: 111 


do 

— + V ■ (pu) = 0 , 
dt 

( continuity ) 

(i) 

dpil _ „ 

-fV*(puu ) + Vp = V- o + pg, 

dt 

( momentum ) 

(2) 

dp h dP 

+ V • (pu h) - — + Q , 

dt dt 

(enthalpy) 

(3) 

P (t)M =p RT , 

w 

(thermal equation of state) 

(4) 

h(T) — f c (i)di , 

D 

(ca loric equation of sta te) 

(5) 


where o is the laminar viscous stress tensor, Q is the volumetric rate of heating 
due to such sources as chemical reaction or divergence of the laminar heat flux. 
Pit) is the volume average pressure of the system, p(x,f) is the pressure fluctua- 
tion from the mean value Pit), and M w is the molecular weight of the gas. For low 
Mach number flows Ipl/P *= M 2 , where M is the Mach number. 1 ? From these 
equations one can derive an equation for the divergence of the velocity field: 1 9 


Vu 


J_ dP 

yP dt 



( 6 ) 


In an open flow system, P is just the ambient pressure; for flow in a closed volume 
V, an equation for P can be derived by integrating (6) over V : 


_L 

yP dt 


1 dV y- 1 J_ 
V dt yP V 


Q dv . 
V 


(7) 


C. The Averaged Equations 

In our turbulence equations we will use both unweighted averaged quanti- 
ties and Favre (density-weighted) averaged quantities. The unweighted average 
and Favre average of a quantity <p are defined respectively by 


- 1 V ' 

<b (x , t) — / <p (x , t) 

v NE^ a 


(8) 


and 
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1 


(9) 


<p(x , t) = - 

p (x , A 


1 


NK 


X P U (X ’ t] ^a U ’ l) 


where ensemble averages are used, <p a being the value of $ in a particular 
experiment a and NE being the number of experiments in the ensemble. The 
fluctuations from these averages are denoted by 


4>' = 4> a - <P 


( 10 ) 


and 


0” = 0 a - 0 , (11) 

where we drop the subscript a on the fluctuations. 

By averaging Eqs. (l)-(5) we obtain the turbulence equations: 


ap - -v 

f V • ( p u) = o , 

dt 

(12) 

ap u — — — 

^ + V • ( p uu) + Vp _ V • R + V * o a- p g , 

(13) 

ap h - — dP — , 

+ V.(p hu)- — + Q - V.^ , 
at at 

(14) 

P M = p RT , 

w r 

(15) 


and 


h = j c(\)d\ , (16) 

where th e Reyn olds stress R is given by - pu’u" and the turbulent heat flux <|>* is 
given by pu”/i\ In deriving Eq. (16) we have assumed that the characteristic tem- 
peratures over which significant changes in c p occur, are much larger than char- 
acteristic temperature fluctuations T 

D. An Alternative Expression for 

An alternative expression for the turbulent heat flux can be derived from 
the averaged and unaveraged equations of state. Subtracting (15) from (4) 
results in 


P'M W = pRT- pRT = pRT" + p'RT . (17) 

We now assume that 

|P'|/ P << ip'l/ p . (18) 
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Under this assumption the left-hand side of (17) can be neglected, and we obtain 


rjitt 



(19) 


wherein temperature and density fluctuations are directly related. 

The assumption (18) requires some discussion. It will certainly be true for 
an open flow system, because the P is always the ambient pressure and never 
fluctuates. In a closed system, such as in internal combustion engine cylinder, 
there can be significant fluctuations in the mean pressure P, primarily due to 
fluctuations in the chemical heat release rate. In an internal combustion engine 
these are referred to as cycle-to-cycle variations, and there is currently some 
debate 20 whether or not these cycle-to-cycle variations should be called turbu- 
lence. We recommend that in performing the averages (8) and (9) one should_use 
only experiments for which the mean pressure history P(t) is nearly equal to P(t). 
When this limited ensemble average is used, the assumption (18) is automat- 
ically satisfied. 

Subtracting (16) from (5) gives 


h" = L c(x)di=c(T)T" , ( 20 ‘ 

] r p p 

where again we use the assumption that c p varies little for temperature changes 
equal to T". Using (20) and (19) in the defining formula for <|>* gives 


** = pu "h" 


= -c ( T)T p'u" 
p 

= p c(T) T a 


( 21 ) 


where 


Equation (21) is the alternative expression we seek for the heat flux. It says 
that <J>^ is proportional to a quantity a that can be loosely thought of as the differ- 
ence between the volume-averaged and mass-averaged velocities. In a two- 
density or two-phase system, a is proportional to the difference between the veloc- 
ities of the two phases. In order to investigate further the nature of the turbulent 
heat flux, we must derive a transport equation for a. 

E. The Transport Equation for a 

The transport equation for a is obtained by subtracting the equation for the 
mass-averaged velocity u from that for the volume-averaged velocity u. The 
result is 
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aa ^ / 1 1 \ - 1 

h u-Vu — u • Vu + ( — — = ) V p + — Vp' 

dt \ p p / p 


= - Iv-R +(- - =)v- o + - V-o' . (23) 

p v P p > P 

As an aid in modeling some of the terms in Eq. (23) we will derive an 
a-equation for two-phase flows and compare this with Eq. (23). It will turn out 
that all of the terms in Eq. (23), with the exception of those associated with the 
fluctuating stresses p' and o', will be d uplicat ed e xactly b y terms in the 
a-equation for two-phase flows. The terms 1/p Vp' and 1/p v • o’ are then associ- 
ated with terms in the two-phase a-equation that arise due to momentum 
exchange between the phases. With this comparison as a guide we postulate a 
model for the fluctuating stress terms. 

We now briefly introduce the equations for two-phase flow. For more details 
the interested reader should consult Refs. 2 and 3. For simplicity we restrict our- 
selves to the flow of two incompressible phases. The continuity equations for each 
phase are 


d Pl Q ! -1 

— +V.(p 1 a 1 u ) = J 21 


(24) 


and 


d P 2 Q 2 

dt 


+ V (p 2 a 2 u) = J 


12 


-J 


21 


(25) 


In these equations p; are the microscopic (or conditional) densities, which we 
assume to be constant; ui are the average velocities within each phase (where the 
"i” next to the overbar indicates a conditional average in phase i)\ a, are the vol- 
ume fractions of each phase (ai + a 2 = 1); and J 21 is the rate of mass transfer per 
unit volume from phase 2 to phase 1. The momentum equations for each phase 
are 


— 1 

6p |C* j u 

dt 


-t-V^PjQjU U 1 + Q^p 


P] a ]g + QjV • o + V ■ (QjRj) + P 2I 


(26) 


and 


d P 2 Q 2 u2 

dt 


- 2-2 


+ V • (p 2 a 2 u u| + a 2 Vp = + a 2 V • a + V • (a^ - P 


21 


(27) 


As is commonly done in two-phase flow modeling, 2 we assume here that the 
phases are in local pressure equilibrium; that is pi = p2 = p. The conditional 
Reynolds stresses R ( are given by 
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H . = -p ( u — u')(u - u ) . 

I l 

The rate of momentum exchange per unit volume from phase 2 to phase 1 is 

denoted by l J 2i- - - : 

Average flow variables f are related to the averages within each phase f by 

f = Q f + a ~T ■ 

Thus, for example, 

p = a , Pl +a 2 p 2 


(28) 

(29) 


and 


-~ -1 -2 
pU=a i P l U + a 2 P 2 U ' 

By using the relations (29) and adding Eqs. (24) and (25) and Eqs. (26) and (27) we 
obtain the same total mass and momentum equations, Eqs. (12) and (13), that we 
have previously derived, when it is realized that 


K = °1 K 1 + °2 R 2 

— 1 ~ “I 

- p d (U - u)(u - u) 

— 2 ~ — 2 

- p 9 a 2 (u - u)(u - u) 

Thus the mass-averaged velocity equation is 


(30) 


du ^ 1 — t — 1 /O M 

— + u • Vu + = = = V ■ a + - V ■ R + g . v<li/ 

dt p p p 

To obtain the a-equation, we will subtract (31) from an equation for u, 
which will now be derived. By dividing (26) by pi and (27) by p2 and summing the 
results, we obtain 


du —1—1 — 2— 2 X / 1 \ — 

— + V • (a.u u + a..u u ) + I - I v p 
dt 1 1 ' p / 


= g + 


(I)v^ + t.(Jk | + Jk 2 ) 


+ P 2l' - 



(32) 


It can be seen that 
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- 1-1 - 2-2 1 2 

a u u + a u u - R - — R„ 

1 1 p 1 o Z 

P 1 f 2 


= V-(uu)= u-Vu + uV-u 


(33) 


Using (33) in (32) and subtracting (31) yields the two-phase a-equation: 


da — — ~ ~ ( \ 1 \ - 

— + u*Vu — u-Vu-h uV-u 4- I - — - V p 

dl ' P p / 


1 

= - = V R + 
P 


1 

P 


= W - o 


p 


+ 1 


. (I 

21 V n 


Comparing (34) and (23), we see they agree if 



(34) 


1 1 / 1 1 \ 

-Vp'- -V -o' = uV • u + P — - — . (35) 

P p 21 V p 2 P, / 

This is the relationship we seek between the fluctuating stress terms and the two- 
phase momentum transfer terms. To obtain closure we need to postulate a form 
for P 21 , and we will investigate expressions employed in two-phase flow 
modeling. 

F. Expressions for the Momentum Exchange Rate in Two-Phase Flows 

For a dispersed phase 2 of equal-sized spherical particles of radius r in a 
continuous phase 1, an expression for the momentum transfer term is21 


P 


21 


3 P 1 Q 2 -2 -1 -2 -1 -2 

— C, |u - u I (u - u ) + J.., u . 

8 u r Zl 


(36) 


Thus, P 21 has two terms — one due to aerodynamic drag and one due to mass 
exchange. This form of the mass exchange term assumes there are no circulation 
velocities within the particles. Equation (36) has theoretical justification l when 
P2 > > pi and velocities within each phase are sharply peaked near their mean 
values. It neglects virtual mass effects, Basset history effects, and particle dis- 
tortions and oscillations. 22 

Motivated by Eq. (36), modelers usually use a similar expression for all two- 
phase regimes: 2 


l* 2l = ^(u 2 - u*) + J, n u‘ s (37) 

where K is called the drag function and u « is some average interface velocity. The 
quantity K is a positive function of pi, p2, aj, 02 , |u 2 - Q1J, and an entity size r. 

If we accept Eq. (37) then one is led to the postulate that the fluctuating 
stress terms in the a-equation (23) contribute to the decay of a. Indeed one can 
show that for a two-phase flow 


a = 


a . a 2 ( P.-P 2 ) -2 

= (u 


p 



(38) 
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and (38), in conjunction with (35) and (37), gives 


1 l ( 1 l 

— V p' — - V - o' = uT u + K 'a + J., I — — — 
P P 2 , \p 2 Pj 

where 


(39) 


p 

K = - K . 

a l a 2 PiP 2 

K has dimensions of a frequency. 

G. Final Form of the a-Equation 

After substituting (39) in (23) and some rearrangement of terms one obtains 

3a ~ ~ / 1 1\ — « - - b — 

- I- V • (ua + au) + K’ a + J I — — — )u +uV-u+=Vd 

91 V P 2 P , 7 P 

1 1 - 1 

= - = V • R — V • (iTu") = = iTu" • V p + - V • ( p'u'u' ) . (40) 

P P P 

Here we have introduced the quantity 6 as a dimensionless measure of the 
density fluctuations: 



If the density fluctuations are not too large, then b is approximately a self- 
correlation coefficient for density fluctuations: 


b = 



In fact. Ref. 4 uses 


(42) 


B = (p') 2 (43) 

as a measure of the density fluctuations. We will develop a transport equation for 
6 in future work. 

Three further terms in (40) must be modeled. First we deal with the mass 
exchange term. One can show from (24) and (25) that 


V- u 




(44) 


and hence the fourth and fifth terms on the left-hand side of (40) combine to give 
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( 45 ) 


uV • u + J ( — — — )u — J ( — — )(u — u) — uV • a . 

2, v P 2 p , ' 21 V p 2 Pj/ 

We assume that U* = u. An assessment of the validity of this assumption must 
await a precise physical int erpr etation of the quantity u s . 

Second, for the tensor u"u" one can show that for two-phase flows 


K, K, 

u"u" = - a. — - a — + 

' P. P 2 

where B is defined by (43) and given in two-phase flows by 



(46) 


B = °i a 2 ( Pi - p/ • 

For future reference we also note that 


(47a) 


6 = — < 47b ) 
PjP-2 

in two-phase flows. We define the volume-averaged conditional Reynolds stress r 

by 


K. «2 

r~a— + a — 

‘ Pi P 2 


(48) 


As a first approximation, and despite experimental evidence to the contrary in 
turbulent flame experiments,^ we assume the conditional Reynolds stresses are 
equal and isotropic. Then 


R. 


R 


2 
P 9 


2 

- - *'l 
3 


(49) 


where k' is related to the specific turbulent kinetic energy k — -y(u*)2 by 


1 a 

k' + 

2 b 


2 


(50) 


A transport equation for k will be developed in future work. 

We also use a two-density distribution to model the triple correlation term 
in Eq. (40). After some algebraic manipulation and use of the assumption (49) 
one obtains 


p'u'u" = - p 



( 51 ) 
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By substituting (45), (46), (49), and (51) in (40) and using the approximation (42) 
we obtain the final form of the a-equation: 


da ~ b - 

— + u- Va + aV-u + a- Vu + V- (aa) + = V p 

dt p 

aa Tp 2 ( ] a 2 \ V p <xo\ 

= — K'a + + - [k - -=?- (52) 

b p 3 V 2 b J p 

H. An Algebraic Closure Approximation 

In numerical computations of multidimensional fluid flows, use of Eq. (52) 
would require solving two or three additional transport equations for components 
of a. Although this is not an unrealistic task for modern computers, considerable 
computational efficiency would result if an accurate algebraic closure approxima- 
tion for a could be found. In this section we present such an approximation based 
on an assumption whose validity must be tested in experimental comparisons. 
The resulting expression for a predicts gradient heat transport, but also contains 
a contribution that predicts the directed flux arising from the interaction of pres- 
sure gradients and density inhomogeneities. 

The assumption we make is analogous to the drift flux approximation of 
two-phase flow modeling. 2 In two-phase modeling this assumption is that the two 
velocity fields are so tightly coupled through the drag terms that characteristic 
drag times are much smaller than characteristic flow times. For us the assump- 
tion is that 


K' >> — , (53) 

L 

where u 0 and L are a characteristic velocity and gradient length for the flow. 

Assuming (53) is true, order of magnitude estimates of the terms in (52) 
show all terms on the left-hand side can be neglected, except the pressure grad- 
ient term. On the right-hand side, the dyadic product term aa is negligible since 
Eqs. (38) and (47) show that a/6 is proportional to the velocity difference between 
fluid elements of different density. The resulting equation for a becomes 


a 


I 

~K' 




Equation (54) can be put in a more recognizable form if we use 


(54) 


Vp 


P 


VT 

~¥ 


a - 


1 

~K' 



(55) 


In conjunction with (21), Eq. (55) gives a heat flux that is the sum of contributions 
proportional to — Vpand — VT. The former is the directed flux. It goes to zero in 
the absence of density fluctuations 6. 
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The gradient transport term in (55) looks similar to the gradient heat flux 
commonly used in turbulence modeling, but there is a difference. The usual form 
used for the turbulent heat flux® is 


<t» = - pc, 


P Pr T t 


. VT 


(56) 


where Pr T is the turbulent Prandtl number and e the turbulence dissipation rate. 
Equation (56) agrees with the heat flux contribution obtained from the second 
term in (55) if the drag time associated with fluid elements of differing density 
equals the turbulence dissipation time. For the momentum exchange function 
(36) it can be seen that the drag time is 


_L ^ P 2 r (57) 

K ' p.iu'-ui 

when cd is approximately unity. On the other hand, the turbulence dissipation 
time is 

I „ L t (58) 

Z k* 

where L is a turbulence length scale. Equations (57) and (58) agree if r ~ L, ~ 

| q 1 _ u2| and pi =» p 2 , but when these equalities are violated more accurate heat 
fluxes could be obtained using a drag time, and not a turbulence dissipation time, 
to evaluate the heat flux vector. 

111. COMPARISON WITH OTHER WORK 

In this section we compare our a-equation with two others in the literature. 
In the BML formulation for turbulent flames, 10 an equation is kept for the tur 
bulent flux of reaction progress variable c. Our quantity a is just a constant times 
the turbulent flux of c: 


pu"c” 


P P 

r p 

P, - P. 


a , 


(59) 


where p r and pp are the reactant and product densities. Two differences are ob- 
served between the a-equation one derives from the BML formulation and ours. 
First, in the BML formulation it is not assumed that the conditional Reynolds 
stresses within each phase are equal and isotropic. An equation for the uncondi- 
tional Reynolds stress R is retained, and the difference between the conditional 
Reynolds stresses is modeled using R. Accordingly, the double and triple correla 
tion terms on the right-hand side of (40) are modeled in a more detailed fashion, 
although the authors observe™ that "this modeling is generally not found to be 
too critical to the predictions of first- and second-moment unconditional 
quantities.” 

The second difference is in the modeling of the fluctuating stress terms. The 
authors follow Launder23 and model 
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(60) 


1 d P' „ e 
- — = 2c, — a — c n 
p dx * c k 2 , 


a • Vu 


where c i c and C 2 c are empirical constants. Comparison with Eq. (39) shows that 
these models have in common the decay of a term and that these would be the 
same if 


K' = 2c, - 
'e k 


(61) 


Besnard, Harlow, and Ra uen zahn 4 keep equations for both the turbulent 
heat flux and the quantity A = p’u', which is related to a by 


since they are interested in more complicated equations of state in which the 
relation ( 19) does not hold. Their a-equation differs from ours in several respects. 
An equation for the Reynolds stress is retained and used in modeling the first 
term on the right-hand side of (40). The triple correlation is broken into two 
terms 


p'u'iT = — 2 p aa + p'u'u' , 

and the latter term is modeled by a gradient diffusion of a. There is a decay of a 
that arises solely from the viscous stress terms. 

IV. SUMMARY AND FUTURE WORK 

We have derived a transport equation for the quantity a, which is the differ 
ence between the volume- and mass-averaged velocities and is simply related to 
the turbulent heat flux Using this equation and an assumption analogous to 
the drift flux approximation of two-phase flow modeling, we have obtained an 
algebraic closure relation for 4>fa that exhibits fluxes due to directed transport 
proportional to — Vp and due to gradient transport proportional to — VT. 

Much work remains to be done before the model can be used in predictive 
calculations of low Mach number flows with large density variations. The equa- 
tion for a involves an additional scalar fa that is a measure of the density fluctua- 
tions. An equation for b must be derived and terms in it modeled. We hope to use 
the a - and fa-equations in conjunction with a k - e turbulence model. The k- and in- 
equations must be reexamined to see what modifications are needed when the 
flows have large density variations. When mass transport is important, such as 
in many combustion problems, expressions for the turbulent mass flux must be 
developed. 

In an effort to test some of the modeling assumptions we are currently 
writing a one-dimensional code that solves the turbulence equations of this 
paper. Computed results will be compared with experimental measurements of 
Rayleigh-Taylor instability, turbulent premixed flames, and flows with centri- 
fuging and density variations. These results and extensions of the model will be 
reported in future publications. 
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